Paper Outline: Snow Interception Processes in the Subalpine
Introduction
There is a need for robust models of snow redistribution by vegetation and wind processes to estimate snow accumulation in mountain forests. To achieve this, a comprehensive understanding of snow redistribution processes is required. However, existing snow interception parameterizations are based on limited observations in distinct climates. Rapid changes in climate and forest ecology illustrate the pressing need to assess whether existing snow interception parameterizations are suitable for predicting snow accumulation in diverse and changing environments.
New observations of snow interception and ablation processes will help determine if existing theories are applicable in mountain forests.
Specifically theories that suggest interception efficiency (interception/snowfall) declines with increasing snow load (Hedstrom & Pomeroy (1998)) or if there is an influence of air temperature (Storck et al., 2002) and wind speed (Hedstrom & Pomeroy, 1998) on interception efficiency will be tested.
The original theory has been simplified over time, i.e. the increase in canopy coverage with increasing wind speed is not included in more recent parameterizations Roth & Nolin (2019).
Updates to theory (e.g., Gelfan et al., 2004) has been ignored in recent studies (e.g., Lundquist et al., 2021).
While many processes occur simultaneously in reality i.e. high winds increasing I/P through horizontal particle trajectory and simultaneously reducing I/P due to unloading. Air temperature may increase I/P through increased cohesion but also may decrease I/P due to melt. Representing these individual processes is important for modelling subcanopy snow accumulation. However since representing all individual processes in models is not ideal it is crucial to determine which are the dominant processes important to represent.
Existing models of snow interception were parameterized using measurements of interception efficiency over snowfall events ranging from hourly (Storck et al., 2002) to weekly timesteps (Hedstrom & Pomeroy, 1998). Based on the different measurement time intervals studies vary in the amount of time possible for ablative processes and may influence model estimates of interception. Despite the inclusion of unloading in the interception parameterizations developed in these studies, they are often combined with additional unloading parameterization in earth system models (Clark et al., 2015; Ellis et al., 2013) leading to some potential of double counting of the ablation process.
Recent observational studies Xiao et al. (2019) do not consider the possibility of subcanopy snowmelt artificially increasing their I/P values.
Uncertainties also arise in the scaling of point or branch scale measurements to the plot scale.
Forest structure governs the interception efficiency observed at a given location (Hedstrom & Pomeroy, 1998; Roth & Nolin, 2019). Metrics used in common snow interception parameterizations (Hedstrom & Pomeroy, 1998; Storck et al., 2002) to describe forest structure include canopy cover and leaf area index (LAI). Leaf area index is defined by Chen et al. (1997) as one half the total green leaf area per unit ground surface area. Canopy cover is defined in Hedstrom & Pomeroy (1998) as the fraction of sky not visible by the instrument from under the canopy. While more detailed forest structure metrics exist derived from detailed LiDAR scans (Helbig et al., 2020; Roth & Nolin, 2019), often they are not available at regional extents required to run hydrological models.
Objective: To examine the dominant processes governing snow accumulation in a subalpine forest, that control the partitioning of snow interception in the canopy and the subsequent unloading to the ground.
Research Questions:
- How is snow interception influenced by meteorology?
- What are the dominant ablation processes observed?
- How does meteorology influence unloading?
Theory
Figure 1 shows the mass exchange of snowfall between the forest canopy and the snowpack on the ground surface. The storage of snow water equivalent (SWE) is represented here with respect to the canopy (\(W\), mm) or the surface snowpack (\(S\), mm). Fluxes that are repeated between the canopy and snowpack control volume have a superscript to specify what control volume they refer to (i.e. \(q^{veg}\) refers to the vegetation control volume and \(q^{snow}\) refers to the surface snowpack control volume).
The change in canopy SWE storage, \(W\) (mm), may be represented as:
\[ \frac{dW}{dt} = q_{sf} - q_{tf}(W) - q_{unld}(W) - q_{drip}(W) - q_{wind}^{veg}(W) - q_{sub}^{veg}(W) \tag{1}\]
where \(q_{sf}\) (mm s-1) is the above canopy snowfall rate, \(q_{tf}\) (mm s-1) is the throughfall rate, \(q_{wind}^{veg}\) (mm s-1) is the wind transport rate in our out of the control volume (typically assumed to be negligible in the literature), \(q_{sub}^{veg}\) (mm s-1) is the intercepted snow sublimation rate, \(q_{unld}\) (mm s-1) is the canopy snow unloading rate and \(q_{drip}\) (mm s-1) is the canopy snow drip rate due to canopy snowmelt. \(q_{wind}^{veg}\) and \(q_{sub}^{veg}\) may be a positive or negative flux. Where all of the above rates are a function of snow load (\(W\)), which is how much snow is present in the canopy.
Estimating Canopy Snow Interception
Current snow interception parameterizations differ in the fraction of snowfall that is intercepted in each storm Figure 2. These parameterizations have all been derived for evergreen needleleaf forests (Hedstrom & Pomeroy, 1998; Moeser et al., 2015; Storck et al., 2002).
These differences are attributed to differing underlying theories, assumptions and processes included within each parameterization. The underlying theory of the Moeser et al. (2015) parameterization stems from the Satterlund & Haupt (1967) study who observed an initial increase in the rate of intercepted snow, as snowflakes bridge gaps between needles. It may also be inferred that during the small near 0°C snowfall events observed in Satterlund & Haupt (1967) the majority of snow may have melted immediately due to a warm canopy. The initial low interception efficiency was followed by an increase and then flattening off of the interception rate as branches bend due to the weight of snow which Satterlund & Haupt (1967) represented by a numerical analytical sigmoidal function. Since the (Moeser et al., 2015; Satterlund & Haupt, 1967) parameterizations were derived on incremental snowfall totals at irregular time intervals, they cannot directly be included in hydrological models. The (Andreadis et al., 2009; Hedstrom & Pomeroy, 1998; Storck et al., 2002) parameterizations were constructed to run on regular time intervals using precipitation rates and are the most common parameterizations used in earth system models (Best et al., 2011; Clark et al., 2020; Ellis et al., 2010; Essery et al., 2016; Verseghy, 2017). These parameterizations will therefore be discussed here in more detail and will be followed by a brief description of the (Moeser et al., 2015; Roth & Nolin, 2019; Satterlund & Haupt, 1967) parameterizations used for snowfall totals.
Hedstrom & Pomeroy (1998)
In the observations by Hedstrom & Pomeroy (1998), snow interception efficiency starts high and then declines and was best represented using an inverse exponential function shown in Figure 2. The Hedstrom & Pomeroy (1998) parameterization therefore differs from the (Moeser et al., 2015; Satterlund & Haupt, 1967) sigmoidal function as it does not include a representation for the initially slow interception rate. However, as shown in Figure 2 all parameterizations include declining the interception efficiency as the canopy fills up with snow.
Using observations in the southern boreal forest from a Jack Pine and Black Spruce stands, Hedstrom & Pomeroy (1998) formulated the conceptual theory to calculate the change in canopy snow storage, \(\Delta W\) (mm) over an accumulation period. In the following equations canopy interception is assumed to occur in isolation from the other processes in Equation 1 and therefore the change in \(W\) over time may be calculated as a modification of Equation 1:
\[ \frac{dW}{dt} = q_{sf}-q_{tf}(W) \tag{2}\]
the parameterization for \(q_{tf}(W)\) is:
\[ q_{tf}(W) = \begin{cases} (1-ip(W))q_{sf},& \text{if } W< W_{max}\\ q_{sf}, & \text{otherwise} \end{cases} \tag{3}\]
where \(ip(W)\) (-), is the interception efficiency which is a function of \(W\). In Hedstrom & Pomeroy (1998), \(ip(W)\) is calculated as:
\[ ip(W) = \frac{\Delta W}{\overline{q_{sf}} \Delta t} \tag{4}\]
where \(\overline{q_{sf}}\) (mm), is the average snowfall rate over the discrete time interval \(\Delta t\). For a given snowfall event, Hedstrom & Pomeroy (1998) calculate \(ip(W)\) as:
\[ ip(W) = C_p \frac{W_{max}-W_0-\Delta W}{W_{max}} \tag{5}\]
the above equation is written in this way in Hedstrom & Pomeroy (1998) since they had measurements of \(W_0\) at the beginning of the storm. However this equation further simplified here since:
\[ W = W_0 + \Delta W \]
and therefore:
\[ ip(W) = C_{lca}\frac{W_{max}-W}{W_{max}} \tag{6}\]
where \(W_0\) (mm), is the canopy snow load before snowfall is added to the canopy, \(\Delta W\) (mm), is the change in canopy snow load due to snowfall, and \(C_{lca}\) (dimensionless) is the snow-leaf contact area ratio. The calculation of \(C_{lca}\) in Hedstrom & Pomeroy (1998) is:
\[ C_{lca} = \frac{C_c}{1-\frac{C_cuH}{wJ}} \tag{7}\]
where \(u\) (m s-1) is the horizontal velocity of the snow particle (approximated by wind speed), \(w\) (m s-1) is the snow particle vertical fall velocity, \(H\) (m) is the height of the canopy, \(J\) (m) is the forested downwind distance.
Hedstrom & Pomeroy (1998) integrate Equation 5 to provide an analytical solution to calculate the change in intercepted snow load, \(\Delta W\) (mm), before unloading which is calculated as:
\[ \Delta W = (W_{max}-W_0)\left(1- exp \left(-\frac{C_{lca}}{W_{max}} \overline{q_{sf}} \Delta t \right) \right) \tag{8}\]
The derivation of Equation 8 from Equation 6 is provided in the appendix.
The maximum intercepted snow load, \(W_{max}\) (mm), is calculated as:
\[ W_{max} = S \cdot LAI \tag{9}\]
where \(LAI\) (dimensionless) is the leaf area index, and \(S\) (mm) is a species maximum snow load correction factor that is a function of snow density:
\[ S = \overline{S} (0.27 + \frac{46}{\rho_s}) \tag{10}\]
where \(\rho_s\) (kg m-3) is the fresh snow density (kg m-3) and Schmidt & Gluns (1991) observed \(\overline{S}\) = 6.6 and 5.9 mm for pine and spruce, respectively.
The calculation for fresh snow density \(\rho_s\) was developed by Hedstrom & Pomeroy (1998) based on observations from sites in Saskatchewan and Yukon (Hedstrom & Pomeroy, 1998), the Fraser Experimental Forest, Colorado (1989) and Nelson, British Columbia (1990) (Schmidt & Gluns, 1991), and the Central Sierra Snow Laboratory, California (Engineers, 1956). Where \(\rho_s\) is calculated as:
\[ \rho_s = 67.92 + 51.25 e ^{(T_{air})/2.59} \tag{11}\]
According to this equation, \(\rho_s\) increases with higher temperatures which will result in lower \(S\) and subsequently \(W_{max}\) Figure 3. In other words, the Hedstrom & Pomeroy (1998) snow interception formula results in slightly lower snow interception efficiency with higher temperatures. The sensitivity of this formula to temperature is seen to be relatively low in the Hedstrom & Pomeroy (1998) analysis. Equation 11 has shown good performance when assessed against 27-years of manual fresh snow density measurements in the Coast Mountains at Whistler Mountain, British Columbia (Barton, 2017).
Storck et al. (2002) and Andreadis et al. (2009)
Andreadis et al. (2009) developed a snow interception model using data collected by Storck et al. (2002) in dense old growth forest in the maritime climate of southwestern Oregon, USA. This method builds off the maximum canopy snow load theory proposed in Hedstrom & Pomeroy (1998) but makes additional modifications to include a step function based on temperature. Here, the snow interception efficiency, \(ip(W)\), from Equation 4 was found equal to a constant of 0.6 based on snow interception observations from Storck et al. (2002) in southern Oregon. Storck et al. (2002) limit, \(W\) as being less than or equal to the maximum interception storage \(W_{max}\) using a step function of temperature.
Here, \(W_{max}\) is given by:
\[ W_{max} = L_r \cdot m \cdot LAI \tag{12}\]
where \(LAI\) is the single-sided leaf area index of the canopy, \(m\) (mm) is an empirical parameter determined based on observations of maximum snow interception capacity (default value = \(5 \cdot e^{-6}\) m), \(L_r\) (dimensionless) is a step function of temperature based on observations by Kobayashi (1987) (on a plank of wood), Pfister & Schneebeli (1999) and Storck et al. (2002) who observed that snow interception decreases with decreasing temperatures. Based on these observations \(L_r\) is calculated as:
\[ L_r = \begin{cases} 4,& \text{if } T_a>-1 ^{\circ} C \\ 1.5 T_a+5.5, & \text{if } -1 ^{\circ} C \geq T_a >-3 ^{\circ} C \\ 1.0, & \text{if } T_a \leq -3 ^{\circ} C \end{cases} \tag{13}\]
The Storck et al. (2002) dataset was derived from observed snowfall events of 5-80 mm with an average of ~15 mm and average air temperature ranging from -10 °C to 5 °C.
A comparison of the Andreadis et al. (2009) and Hedstrom & Pomeroy (1998) canopy interception canopy interception capacity parameterizations is shown in Figure 3.
Methods / Study Area
In-situ observations of forest snow accumulation, snow interception, and canopy snow ablation were collected from at Fortress Mountain Research Basin (FMRB), AB ~ 51 ◦ N, 2100 m asl., a continental headwater basin with sparse subalpine forest (Figure 4 (a)).
A flux tower provided measurements of air temperature, relative humidity and wind speed at 15 m and 7 m above the ground. A weighing precipitation gauge (Pluvio) provided measurements of the snowfall rate (\(q_{sf}\)) to an open clearing adjacent to the study area shown in Figure 4 (b). Throughfall was measured over a discrete time interval (\(\Delta {tf}\)) using areal fresh snow surveys. Throughfall was calculated as the difference of pre and post snowfall event measurements or fresh snow accumulated above a snowpack crust layer with a known date. Continuous fifteen minute measurements of throughfall was also measured using subcanopy-lysimeters shown in Figure 5 (a) and locations shown in Figure 4 (b). During events where the ablative processes in Equation 1 could be considered negligible, Equation 2 was used along with measurements of snowfall and throughfall to infer the amount of snow intercepted in the canopy (\(W\)). A weighed tree lysimeter shown in Figure 5 (b) was used to calculate the weight of snow intercepted in the canopy. The weight of snow in the canopy was scaled to an areal estimate of canopy snow load (\(W\), kg m-2) using the manual snow survey measurements (as in Hedstrom & Pomeroy (1998)).
Results
Forest Structure
During events with lower wind speeds I/P is strongly associated with forest structure metrics LAI and CC across all of the snow survey stations. Canopy closure explains the most variability in I/P at a zenith angle of 30 degrees, while LAI has more variability explained at 35 and 60 degrees. A wider zenith angle incorporates more information about the forest surrounding the snow survey station and may better represent the forest structure important in controlling snowpack accumulation in mountain forests with high wind speeds which lead to horizontal snowflake trajectories. During high wind speed events (10 m wind speed > 3.9 m s-1) the explanatory power of LAI and CC is reduced and an increase in I/P is observed at locations with low LAI and CC while reducing I/P in locations with higher LAI and CC due to wind induced unloading. Since snow survey data was used here some unloading and wind redistribution may have also affected these results.
Meteorology
Air temperature was observed to have little influence over interception efficiency measured using the weighed tree and lysimeters. Evidence of a positive influence of warm temperatures raising I/P was observed with the snow survey data however after careful investigation of time-lapse imagery this increase was due to subcanopy snow melt and not due to snow caught in the canopy which artificially raised the I/P value.
Interception efficiency increased between wind speeds of 0 and 2 m s-1 due to an associated increase in canopy contact area, while later decreasing intercepted load due to increased snow unloading. The increase in I/P with increasing wind speed (below 2 m s-1) is expected as the horizontal particle trajectory of the snowflake becomes more horizontal thereby increasing the effective canopy coverage of the forest.
Interception efficiency increased as the canopy filled with snow and declined later at higher loads.
Canopy Storage Capacity
Existing observations have suggested a relationship in the peak canopy storage capacity and air temperature. Figure 8 shows the increase in canopy storage prior to unloading for 25 snowfall events across a range of average event air temperatures ranging from -20 to near 0 °C. The rate of increase in canopy storage or the maximum canopy storage reached did not exhibit a systematic trend with average event air temperature. No significant relationship was found between the maximum canopy storage reached within event before unloading was plotted aginst average event air temperature (Figure 9).
Discussion
With increased wind speed there is a corresponding increase in canopy contact area due to falling snow particle trajectories in the wind becoming more horizontal. During snowfall events with higher wind speed the forest strcutre metric for a snow survey station was less related to the observed interception efficiency at that location. This is due to the surrounding forest structure intercepting snowfall as the trajectory becomes more horizontal. This thoery is included in Hp98 but is often neglected in models, since it has yet to be confirmed with observed data.
Existing theory suggests either a positive (Andreadis et al., 2009) or negative (Hedstrom & Pomeroy, 1998) of air temperature and the maximum intercepted snow load (Figure 3). The negative association between air temperature and \(W_{max}\) in Hedstrom & Pomeroy (1998) is based on week long snow survey measurements and likley included some ablation (sublimation + unloading). The positive association between air temperature and \(W_{max}\) in Storck et al. (2002) may be a result of the limited frequency of observations at colder temperatures. The new observations collected from FMRB suggest there is no relationship between air temperature and maximum canopy storage (Figure 9) or interception efficiency (Figure 7).
- Some evidence for initial increase in I/P with increasing canopy snow load (Moeser et al., 2015; Satterlund & Haupt, 1967).
Conclusion
- Forest structure metrics calculated using high zenith angles better described the variability in interception efficiency.
- Interception efficiency increased with increasing canopy coverage however the strength of this association was reduced at higher wind speeds.
- High wind speeds increased interception efficiency due to an associated increase in canopy contact area, while later decreasing intercepted load due to increased snow unloading.
- No influence of air temperature on interception efficiency was observed.
- Interception efficiency increased as the canopy filled with snow and declined later at higher loads.
Appendix
The derivation of Equation 8 from Equation 6 is provided by first combining Equation 4 and Equation 6.
\[\frac{1}{W_{max}-W} dW = \frac{C_{lca}}{W_{max}} q_{sf} dt\]
\[\int^{W_1}_{W_0}\left[\frac{1}{W_{max}-W}\right] dW =\frac{C_{lca}}{W_{max}} \overline{q_{sf}} \int^{t_1}_{t_0} dt \]
\[-ln\left(\frac{W_{max}-W_1}{W_{max}-W_0}\right) = \frac{C_{lca}}{W_{max}}\overline{q_{sf}} \Delta t \]
\[\frac{W_{max}-W_1}{W_{max}-W_0} = exp \left(-\frac{C_{lca}}{W_{max}}\overline{q_{sf}} \Delta t \right)\]
\[1-\frac{W_{max}-W_1}{W_{max}-W_0} = 1- exp \left(-\frac{C_{lca}}{W_{max}}\overline{q_{sf}} \Delta t \right)\]
\[(W_{max}-W_0)-(W_{max}-W_1) = (W_{max}-W_0)\left(1- exp\left(-\frac{C_{lca}}{W_{max}} \overline{q_{sf}} \Delta t \right) \right)\]
\[\Delta W = (W_{max}-W_0)\left(1- exp \left(-\frac{C_{lca}}{W_{max}}\overline{q_{sf}} \Delta t \right)\right) \tag{14}\]
here, it is assumed that \(\overline{q_{sf}}\) is the average snowfall rate over the discrete time interval \(\Delta t\). Since \(\overline{q_{sf}}\), \(C_{lca}\) and \(W_{max}\) are temporally constant over the discrete time interval \(\Delta t\) they can be moved outside the integral. The analytical solution in Equation 8 is only possible because canopy snow interception is treated in isolation from the other processes in Equation 1.
The exact analytical solution in Equation 8 is a typical method used in hydrological models to discretize parameterizations (e.g. Equation 6) over discrete time intervals. However, problems with this method have been outlined in Clark & Kavetski (2010) and Schoups et al. (2010) and attribute approximation in time stepping schemes as the source of considerable model error.